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Abstract 

A new class of models based on hysteresis functions is developed to de- 
scribe atomic force microscopes operating in dynamic mode. Such models 



are able to account for dissipative phenomena in the tip-sample interac- 
tion which are peculiar of this operation mode. The model analysis, which 
can be pursued using frequency domain techniques, provides a clear in- 
sight of specific nonlinear behaviours. Experiments show good agreement 
with the identified models. 

1 Introduction 

Physical systems showing impact phenomena are frequent in many fields [T]. 
Main applications occur in mechanics where macroscopic objects are consid- 
ered. In this case, an impulsive approximation for interaction forces with a pure 
repulsive nature can often be correctly assumed. Moreover, energy losses are 
traditionally considered by introducing the concept of coefficient of restitution 
[2] [3]. However, there are many situations when this kind of approximation 
can not be considered satisfactory, for example when the interaction involves 
both attractive and repulsive forces or when the interaction can not be assumed 
instantaneous. The aim of this work is to exploit a hysteresis function to model 
the related interaction forces. This new model can be viewed as a generalization 
of the impulsive case and allows for the use of potential functions even if the 
system is dissipative. It also presents advantages when the interaction forces 
involve both repulsive and attractive parts or when the duration of the impact 
is not neglegible. In addition the hysteresis model allows for the use of powerful 
analysis techniques, such as harmonic balance [4], which could not be used for 
impulsive forces. To show how the above impact model can be successfully em- 
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ployed, its application to an Atomic Force Microscope (AFM) is demonstrated. 
Specifically, we limit ourselves to the study of an AFM operating in dynamic 
mode, whose schematic is depicted in Figure [T] the cantilever is periodically 
forced by a piezo placed under its support inducing a periodic oscillation that 
is influenced by the interaction forces between the cantilever tip and the sam- 
ple. The topography can be inferred by slowly moving the cantilever along the 
sample surface by means of a piezoactuator and by measuring the amplitude of 
the cantilever deflection through an optical lever method. A feedback controller 
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Figure 1: Schematic of a tapping- mode AFM. 



driving the piezo input voltage is employed to reject variations of the separation 
between the sample and the tip center of oscillation due to the sample topog- 
raphy. In AFM, the amplitude-distance curve is not used to obtain the sample 
topography. It is the control-signal that is used to get the image. Therefore, the 
amplitude-distance curve for topography is not crucial. However, topography is 
not the only information one might be interested in. One of the primary uses 
of AFM is the study of force interactions [5j. Two methods are prevalent. The 
cantilever-sample offset (also termed as separation), which is a measure of the 
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distance of the cantilever holder and the piezo-positioncr, is first reduced by 
the piezo positioner, where the sample surfaces approaches the cantilever-tip. 
In the retract phase the cantilever-sample offset is increased by using the piezo 
positioner. During the approach and the retract phases the cantilever deflec- 
tion signal is recorded. The force felt by the cantilever, can be obtained by 
multiplying the deflection by its spring-constant. By plotting the force felt by 
the cantilever against the cantilever-sample offset the force curves are obtained. 
These curves are called static force curves. In dynamic force curves the can- 
tilever is oscillated using the dither piezo. The amplitude of the first harmonic 
is plotted against the cantilever-sample offset during the approach and retract 
phases. The dynamic force curves are gentler on the sample and therefore are 
the preferred means of investigating samples that are soft (e. g. biological 
samples). One of the difficulties of using the dynamic force curve mode when 
compared to static force curve mode in obtaining force curves is that determin- 
ing force-separation curves from the measured amplitude-separation curves is 
not as straightforward. In most cases, dynamic force curves are obtained by in- 
tensive numerical simulation. For example, in |B] and [7] models that accurately 
describe the device behaviour are proposed. In another approach, parametrized 
models of the tip sample interaction are assumed, the parameters identified us- 
ing the amplitude-separation data, and subsequently, the force-curve data is 
generated using the identified model. In 8 an identification algorithm of the 
force-curve is obtained by the numerical computation of an explicit integral 
equation. There are only few attempts in the literature on analytical results. 
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Analytical results can be found in [S^, where a simple impulsive impact model 
is developed. However, since the employed model neglects attractive forces, it 
does not seem able to explain some important characteristics of the tip-sample 
interaction observed in experiments. 

In this work, we develop a complete frequency analysis of a dynamic-mode AFM 
exploiting the proposed hysteresis description and taking into account attractive 
forces in the sample-cantilever interaction. The main feature of the proposed 
model is to provide results without the means of numerical simulations, for 
example evaluating the separation-amplitude curve for a large class of interac- 
tion forces comprising some of the common potential functions studied in the 
literature, such as the classical Lennard- Jones potential lOj. Other peculiar 
attractive features of the proposed class of models are: i) it can easily account 
for energy losses; ii) it is suited for nonlinear frequency-domain identification 
techniques such as those proposed in 7J and [11 ; in) it facilitates to study some 
structural properties of the system such as bifurcation phenomena experimen- 
tally observed exploiting frequency domain techniques as in [12j . Identification 
results based on experimental data are provided where the hysteresis model 
gives a good qualitative and quantitative characterization of the tip-sample be- 
haviour. 

The paper is organized as follows. In Section [2] we briefly describe the gen- 
eral problem of modelling an impact reminding many consolidate notions for 
the sake of clarity. In Section [3] we exploit such a model to describe the AFM 
tapping-mode dynamics and in Section [4] a frequency analysis is provided us- 
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ing harmonic balance techniques. In Section [5] the identification procedure is 
described and finally in Section [S] experimental results are discussed. 



2 Hysteresis functions to model a collision 

Let Pi and P2 be two material objects with masses toi and m2, respectively, 
moving along the x axis, with position xi and X2 {xi < X2)- We consider Pi 
and P2 subject to external forces fi{t) and /2(i) respectively, and to a mutual 
internal force. Pi exerts a force on P2 given by /12 and P2 exerts an equal 
and opposite force hi on Pi. The interaction force hi are dependent on time 
t, relative separation xi — X2 and relative velocities xi — X2. The following 
dynamical relations describe the system 



Earlier interaction models usually neglected dissipation losses or used a constant 
coefiicient of restitution to account for such losses. Defining 6 :^ X2 — Xi, we 
suppose the interaction between the two masses is negligible outside a time 
interval [is,^/] where S{ts) = S{tf). We intend to limit our study to the case 
where Pi and P2 get closer at the beginning, and then further. We assume that 
the system dynamics can be split in two different phases: an "approach phase" 
in the time interval [tg , f| where S < and a "retract phase" in the time interval 
[t, tf] where 5 >0. We also consider that 6{t) is a continuous function and that 
the set of points where 6 — has zero measure. The model we intend to employ 
in this work defines a particular form for the interaction forces and, at the same 



miXi = fi{t) + hi{t,X2 ~ Xi,X2 - Xi) 



(1) 



m2X2 = f2{t) + h2{t, X2 - Xi,X2 - Xl)- 
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time, allows one to generalize the case of constant coefficient of restitution, not 
only for an instantaneous impact time. In addition, it presents advantages in 
the study of impacting systems with periodic behaviours. The interaction force 
assumes two different forms during the approach and the retract phases 



h{6,6) := < 



h+(S) ifj>0 

(2) 



h-{S) i{S<0. 

The simplicity of the dependence of h on the sign of 6 leads to a tractable 
analysis while capturing the prominent features of finite time that can have 
both attractive as well as repulsive forces. 

If /i+ and h~ are integrable, then the potential functions and U~ can be 
introduced with 



As h-{S) >h+{6) V (5 < S{ts), we have also that U-{S) >U+{6) V 5 < 6{ts). 
Considering the instant i when the relative distance S{t) is the smallest, we can 
state that the potential interaction energy is U~{t) if t <t, while it is U~^{t) if 
t > i. At any instant the interaction force is conservative except at t = i when 
the relative velocity is zero and we have an "instantaneous" energy variation 
AE equal to 

AE = U-{5{t))-U+{S{t)). (4) 

The energy lost in the impact can be interpreted as the area between the curves 
h~{S) and h'^{S) in the interval [5{t),5{ts)]. In the next section we will show 
how this class of hysteresis functions can be exploited to analyze and identify 
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tip-sample interactions in AFMs. 



3 AFM model 

AFM cantilevers can be modeled as a feedback interconnection of a linear system 
L and a nonlinear static function h as depicted in Figure [21 Models with this 
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Figure 2: A feedback interconnection of a linear system and a nonlinear static 
function. 

peculiar structure are well-known as Lur'e models [4]. The system equation can 
be conveniently written using the symbolic form 

yW = i(|)[M'5W,^W) + 7W] (5) 

where y(t) is the measured output (that is the cantilever tip deflection), apart 
an additive constant, is the separation. Sit) = yit) + 1 represents the tip-sample 
distance, and 7(f) is the external periodic forcing 

7(t) = rcos(wt + (j)) . 



8 



The subsystem £ describes the free cantilever dynamics, whose frequency re- 
sponse L{iuj) can be precisely identified using thermal noise or a simple fre- 
quency sweep excitation when the sample is absent |13j . 

The feedback subsystem h accounts for the sample interaction force, which is 
a highly nonlinear function of the tip-sample distance S. Modeling h is still a 
challenging task. The main difficulty lies in the choice of a suitable class of 
functions to describe the force potential. It is a common choice to consider h as 
the sum of a conservative force hcon and a dissipative one hdis 

h{6,5)^Kon{S) + hd^s{.S/5) (6) 

giving to hdis a simple form to allow easy computation P^. Also in [Tl^ it is 
proposed 

Kus{5,5)^V{5)5 (7) 
where F represents a sort of damping coefficient. 

In this paper we consider the following class of hysteresis functions which gen- 
eralizes the one presented in |15| 

E K~K{5) if ^ < 

(8) 

E KtK{5) if > 

n=l 

where hn{S) are a class of suitable non- negative functions where the depen- 
dence on S occurs in h considering only its sign as described in the previous 
section. Relation ([5]) represents a vector space of hysteresis functions made of 
two different positional forces: the first one acts when the tip and the sample 
are approaching and the second one when they are getting further. In order to 
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make the interaction described by h(S, S) dissipative, some constraints on the 
parameters and can be imposed. For example, the condition 

K- > K+ (9) 

makes every base element dissipative. 

In our analysis, we will consider two special cases of this hysteretic interaction. 
This way of modeling dissipations has already been proposed in [16j and has 
been successfully exploited by [15| in an identification procedure. 

3.1 Piecewise interaction force 

The first class of potential functions we treat contains the functions h{-) in the 
form ^ where TV = 2 and 

{0 if (5 > 

(10) 
if ,5 < 0. 

In Figure [3] the shape of such a kind of functions is depicted. Since 6 = y + I, 
here the parameter I models the cantilever deflection at which the tip-sample 
interaction forces become effective. 

3.2 Lennard- Jones-like interaction force 

The Lennard- Jones potential 

f'i^) = ^ + na < rv e AA; X„„ , K^^ e TZ (11) 

is a common choice when fitting statically measured curves often used as a model 
of interaction potential between atoms (see [5j) We consider a generalization of 
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Figure 3: Interaction force modeled by a piecewise linear function. 

the Lennard-Jones Potential in the form ([8]) where 

1 



hn{y) 



(5" ■ 



(12) 



The parameter I represents the cantilever deflection where the interaction force 
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Figure 4: Sketch of a Lennard- Jones-like interaction with hysteresis. 



becomes infinitely large (Figure H]). The choice of this class of functions is 



11 



motivated by its semplicity and also by the fact that long-range dissipative 
interfacial forces has already been successfully modelled using a time-dependent 
power law where the strength of the force depends only on whether the probe 
approaches or retracts away from the sample. [15] 



4 Frequency Analysis via Harmonic Balance 

The linear part of the Lur'e system given by £ in ([51) typically shows a sharp 
filtering effect beyond the first resonance peak because of a very high quality 
factor of the cantilever. Indeed, it is experimentally observed that the cantilever 
trajectory has a quasi-sinusoidal behaviour. The cantilever-tip motion can be 
approximated by 

y{t) ~ yi (t) ■.= Re[A + Be*"'] = A + B cos(wi) (13) 

The corresponding output of the nonlinear hysteresis block can be approximated 
as 

hiy + l,y)^ h{yi + I, yi) ~ Re [NoA + TViBe*"*] (14) 

where 



1 1 '•^ 



No = No{A,B,Lo):=j- / h{yi{t) + l,yi{t))dt 





1 2 '•^ 



m=NiiA,B,Lj) := -- / h{yiit),y,it))e~'^'dt 



(15) 







are the constant and harmonic gains of the nonlinear block also known as the 
describing functions of the nonlinearity 4J. We remark that A^o^ and NiB 
are the first two Fourier coefficients of h{yi{t) + l,yi{t)), thus expression (|15p 
represents a first order harmonic truncation. For the general class of hysteretic 
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force models introduced, we obtain 



No 
" 2 



— V / K+K{1 + A + Bcost)At+ K-K{1 + A + B cos t)At 



hn[B{q + C0ST)]dT 



2 

where 
and 



Similarly, we find for iVi 



E„ := K- + K+ (16) 



A^i = — ^ f / K+K[B{q + cosT)]e-*"dr + / K~K[B{q + cosT)]e-'"dr 

n— 1 ^ ^ 

1 / /■+'^ r^"^ 

= ( S„ / /i„[B(g + cos r)] cos rdr — iA„ / /i„[i?(q + cost)] sin rdr 

T^B V Jo Jo 



(18) 

where 

A„ := K+ - K- . (19) 
Substituting in ([5]), assuming a sinusoidal forcing ^{t) — Re[re*'^'^*+'^)], yields 

A + Be*'^* = -L(0) A^o^ + i(iw) [-N^B + re^-^Je'"* (20) 

or, equivalently, 

[l + L(0)A^o(A S)]A = 



(21) 



[1 + L{iujt)Ni{A, B)] B = L{iujt)Te 



.i<t> 
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Finally, we can easily decouple the variable <j> from (PT|) as follows 



[l + L{0)No{A, B)]A = 

\l + L{iujt)Ni{A,B)\B ^ \L{iujt)\r. (22) 
(/) = arg [LiiLu)-^ + NiiA, B)] . 
The equations (I22p represent a system of three nonlinear equations in the three 
unknown A, B, (p. By solving it, we can find the sinusoidal approximation of 
y{t) given by |T3]). 



4.1 Piecewise interaction model analysis 

For the piecewise- linear potential described in Section [5TT] we obtain 
iVo- i[Eii?i(q) + E2i?2(g)S] 

where 



acos((7) qa.cos{q) - y^l - 
^nQ) — n2(q) := 



2tt 



27r 



^1(9) 



TT 



TT 



S2iq) 
T2{q) 



acos(g) - g^l - 
27r 

27r 



Finally, by the substitutions 



:=Sii?i(<z) r!(g) := S2i?2(<z) 
:= ^iSi{q) + iAiTiiq) £3^2(9) + iA2r2(q) 

we can obtain for the describing functions 

No^j[x{q) + mB] 



(23) 
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In this model, the variable q represents the "penetration" of the tip into the 
sample. In fact, assuming as exact the first harmonic approximation, we have 
that for q > 1 the tip does not get in contact with the sample; for q — I the 
tip grazes the sample and for q < 1 the tip enters the sample. The case q < —I 
does not have a physical meaning in this model. From (PT|). it is also possible 
to write i? as a function of q. In fact 

r 

^ " |LM-i+iVi| ^^^^ 

implies 

iLiiioy^B + BNif = \L{iuj)~''^B + <i> + <fBf = (25) 
The substitutions l> (f> and '5 := * + L{iuj)~^ yield 



l^pB^ + 2Re[$**]B + - r2 = 0. 
which is a simple second order algebraic equation whose roots are 



(26) 



-Ref^^-*] ± jRe[$^'*F - |*|2(|$|2 _ r2) 
B{q) = (27) 

Substituting in (|22|) and reminding that I — qB — A, we can finally write 

A{q) = -L(0) [x{q) + ^{q)B{q)] 

m = arg [L{ioj)-^ + NM. B{q))\ (28) 
l{q) = qB{q) + i(0) [x{q) + ^{q)B{q)] . 
The variable q P?|) . depends on A, B and I (see ([IT])), therefore equations 
(P5|) are implicit relations. System ([^5]) can not be solved in closed form since it 
involves transcendental equations. However, it is possible to obtain its solution 



15 



through a conceptually easy method. Assuming that I is a known parameter 
of the model, it is possible by the last of (|28p . to determine the corresponding 
values of q and then A, B and by exploiting the remaining equations. In other 
words, we have transformed the problem of solving the whole system (1211) into 
the easier problem of solving a single real equation in the unknown q. 
Experimentally, the separation-amplitude curve is obtained by slowly moving 
the sample towards the cantilever and measuring both the amplitude of the first 
harmonic and the separation. Although it is not possible to derive an explicit 
analytical form for B = B{1), we can give a parametric form for it. By using 
the "g-explicit" equations in (|28p we can consider the parametric curve 



4.2 Lennard Jones-like hysteretic model analysis 

For the generic hysteretic interaction force of the class (|12|) . we can evaluate the 
describing functions iVo and A^i of the nonlinearity h: 



I = l{q) 



< 



V g G R. 



(29) 



B = B{q) 




where the functions 






can be analytically evaluated for any given n and q > 1. 
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Imposing harmonic balance, we get 




(30) 



The second equation of ([50]) can be expressed in the form 



N 



L{ioj)Te"l' = L{iuj) 



+ B. 



(31) 



We can remove (j) by multiplying each term by its conjugate. Finally, multiplying 
by B^^ the equation can be easily rewritten as a {N + 2)-degree polynomial in 
the variable B whose coefficients depend only on the variable q 



It can be shown that C2N+2 = C2N+1 = and C2N — — T^. For 

sufficiently large q (that is when the interaction is negligible) we have that 
Ck ^0,yk< 2N, therefore 



One root of the equation above is i? = T\L{iijj)\. This solution corresponds to 
the free oscillation amplitude that the cantilever assumes when the sample is 
far away and does not influence the cantilever dynamics. For every q > I, the 
polynomial equation (|32p can be solved in B. Only the solutions that are real 
and positive have relevance. The constant component of the periodic solution A 
can be evaluated exploiting the first of pO|) . The phase can also be similarly 
obtained as a function of the parameter q 



2N+2 



p{B) = Cn{q)B^ = 0. 



(32) 



p{B) ^ {\L{iuj)\-^B^ - T^)B^^ = 0. 



(33) 




(34) 



17 



Finally, the parameter / is given by the original relation 

l{q)^qB{q)-A{q). (35) 

The final result is that the variables A, B, (p and I are all expressed with respect 
to the parameter q. The separation-amplitude diagram can be obtained con- 
sidering the pair {l{q), B{q)) which describes a curve in a parametric form. A 
similar procedure can be used to obtain the relation between any two variables 
with no need of simulation tools. 

In [T7] and [TB] it is shown that the approximation error of the HB method for 
the analysis of this model is negligible when compared to results obtained by 
simulating the same model. 

5 Identification of the Tip-Sample Force Model 



In this section we present methods to identify parameters of the hysteresis based 
models developed in earlier. As previously discussed, the frequency response 
L{iuj) is known since it can be independently estimated. We also assume that 
the separation I can be changed by means of the piezoactuator placed beneath 
the sample. Therefore, we can consider a set of M experiments with different 
values of I 

Im ■= h + iTT-d rn = l,...,M (36) 
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where Iq is a fixed offset and d > is a suitable separation step. For every 1^ 
the quantities Am, Bm, 4>in can be evaluated from the measured signal y{t) after 
it has reached its steady state, and qm can be computed from (jl7p . 
The functions ([8]) chosen to model the interaction have the useful property that 
it is linear in the parameters K~ and K^, or in virtue of 1161 and I19[ S„ and 
A„. The linear dependence on the parameters aids their identification using the 
harmonic balance relations (|22[) . The first order harmonic balance equations 
lead to a set of M linear equations in the 2N unknown variables I]„ and A„ 



rcos(0,„) - Im[L-i(fc^)]B,„ = S„%^ 
rsin(,/)„0 ~ Re[L-\iLo)]Bm = En=i 



m=l,...,M (37) 



Assuming that there are M > 2N experimental measures and adopting a more 
compact notation, we can write two independent matrix equations 



Ps^ = Qs 
PdA = Qd 



where 



S := 



V J 



A := 



are the unknown vectors and 



(38) 



(39) 



Ps[m,n] Qs[m] := rcos(0,„) - Re[L-\iLj)]Br, 



PD[m,n] := -^iSi^ Q D[m] := T sm{(l),n) - lm[L ^{iu;)]B„ 



(40) 



are constant matrices. 

Since the number of equations is greater than the number of unknowns, (|38[) is 
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not expected to be feasible. A common strategy is to find the set of parameters 
which better "fits" the equations according to the quadratic cost function 

FGo, S, A) = WQs - PsSf + WQd - PoAf, (41) 

where we have stressed the dependence on the offset since it is not apriori 
known. Thus, the optimal values S and A can be evaluated casting an opti- 
mization problem which also takes into account the constraints ^ 

(SGo),A(Zo)) = argminyao,5],A) 
subject to 

(42) 

A < 0. 

We remark that forcing the condition A = in (|42)) is equivalent to the assump- 
tion of a interaction force with no hysteresis and therefore conservative. 
Given Iq, problem (|42p is a quadratic optimization problem with linear con- 
straints. Many algorithms are known in literature to determine its solution 
(S(/o), A(Zo)) ^18j. Finally, we can estimate the offset Iq by solution of the 
following problem 

^0 = argminl^(?o,S(/o), A(?o),) (43) 
which is another minimization over a scalar variable, solvable using a grid strat- 
egy. The identified parameters are (S(Zo), A(Zo)). 

6 Experimental Results 

An atomic force microscope was operated in dynamic mode using a silicon can- 
tilever of 225 fim in length. Using a thermal-response based approach the can- 
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tilever has been identified with a second order linear oscillator with natural 
frequency tOn = 27r73.881 rad/s and damping factor ^ = 0.00378. The spring 
constant K of the cantilever has been estimated to be AnN /nra. A sinusoidal 
voltage with frequency aj„ has been applied to the dither piezo in order to make 
the cantilever oscillate. Experiments were performed on a silicon wafer sample. 
The separation amplitude curve has been experimentally measured during the 
two phases: approach and retract. Both curves show a jump phenomenon oc- 
curring at two close but different values of the separation (dashed and dotted 
curves in Figure [51). Such phenomena are present and documented in literature 
[B] . The identification technique described in the previous section has been em- 
ployed using the data obtained during the retract phase only, while the data 
acquired during the approaching phase have been used for validation purposes. 
The results obtained using the piecewise linear interaction model are reported 
in Figure O (solid curve). The retract curve is well explained by the model data; 
this is not surprising as the data used to obtain the model parameters is the 
retract phase data. For the approach curve, a jump phenomenon occuring at a 
different separation is qualitatively well-predicted, but it can be argued that it 
is not quantitatively satisfactory. 

As a second case, the following simplified model than ([8]) is used 



that is a standard 6—12 Lennard-Jones potential function with a hysteresis 
dissipation. Identification results are shown in Figure [5] As it is evident from 




if (5 > 



if (5 < 



(44) 
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the figure, the model predicts the discontinuity in the approach phase of the 
force curve accurately. The phase of the first harmonic is also predicted well by 
the model. 

As remarked in Figures [5] and [6l harmonic balance has also allowed to reveal 
the presence of instable periodic orbits (in the region in between the two jump 
points) and to clearly explain bifurcation phenomena in the system. 

7 Conclusions 

In the paper we have proposed a class of models for tip-sample interaction in 
atomic force microscopy via impact dynamics. The use of a hysteresis can be 
well combined with harmonic balance techniques for the analysis of oscillatory 
behaviour to provide interesting insights into the dynamics. For instance, the 
presence of jump phenomena discovered in many experiments is well-predicted 
and explained. The suggested method is based on a first order harmonic ap- 
proximation and gives good quantitative results since the linear part of the con- 
sidered Lur'e system shows a sharp filtering effect near the resonance frequency. 
In such a situation, the Harmonic Balance technique has advantages over stan- 
dard numerical approaches since it requires a computational effort much smaller 
than the one required by simulation tools. 
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Figure 5: Experimental separation-amplitude (a) and separation-phase (b) 
curves fitted using the piecewise linear model for the interaction force. Solid 
line is the curve obtained by the model; the dashed and the dotted ones are the 
experimental approach and retract curves respectively. 
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Figure 6: Experimental separation-amplitude (a) and separation-phase (b) 
curves fitted using the Lennard- Jones model for the interaction force. Solid 
line is the curve obtained by the model; the dashed and the dotted ones are the 
experimental approach and retract curves respectively. 
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Figure 7: Identified dynamic Lennard- Jones force curve with hysteresis. 
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